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Abstract 

On the basis of gravity and radar observations with the Cassini spacecraft, the moment of inertia 
of Titan and the orientation of Titan's rotation axis have been estimated in recent studies. According 
to the observed orientation, Titan is close to the Cassini state. However, the observed obliquity is 
inconsistent with the estimate of the moment of inertia for an entirely solid Titan occupying the 
Cassini state. We propose a new Cassini state model for Titan in which we assume the presence 
of a liquid water ocean beneath an ice shell and consider the gravitational and pressure torques 
arising between the different layers of the satellite. With the new model, we find a closer agreement 
between the moment of inertia and the rotation state than for the solid case, strengthening the 
possibility that Titan has a subsurface ocean. 

1 Introduction 

On the basis of Cassini radar images, [6] and [7] precisely measured the orientation of the 
rotation axis of Titan. Using the orientation of the normal to the orbit of Titan given in the IAU 
recommendations (Seidelmann et al. 2007), they determined the obliquity to be about 0.3°. They 
also showed that the rotation axis makes a small angle of about 0.1° with respect to the plane 
defined by the normal to the Laplace plane and the normal to the orbit and, as a result, they stated 
that Titan is close to the Cassini state. 

From Cassini radio tracking, the quadrupole field of Titan has been found to be consistent with 
a body in hydrostatic equilibrium with a moment of inertia C '/MR 2 = 0.3414 ±0.0005 |4j- However, 
in a study of the Cassini state generalized to a multi-frequency orbital node precession in which 
Titan was considered as an entirely solid body, [2] found that the 0.3° obliquity implies a moment 
of inertia of C = 0.55A1R 2 , which is not only inconsistent with the above hydrostatic value of the 
moment of inertia but would also imply a physically implausible interior structure with higher mass 
density towards the surface than towards the center. In their conclusion, the authors proposed that 
a liquid ocean could partially decouple the shell from the interior. 

The obliquity of the Cassini sate, which is an equilibrium orientation of the rotation axis of 
a synchronous satellite, can be derived from angular momentum equations in much the same way 
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Table 1 Orbital theory of Titan and solid Cassini state : columns 2 to 5 : Amplitudes, 
frequencies, periods, and phases of the orbit precession adapted from [TT]. They give the orbital 
precession in the equatorial plane of Saturn and used J 1980 as the time origin. Here we consider the 
Laplace plane and J2000 time origin. The Laplace plane has a node of 184.578° and an inclination 
(also called tilt) of 0.6420° with respect to the equatorial plane of Saturn. Since the tilt is small, the 
frequency and amplitude of the orbital precession are almost the same with respect to the Laplace 
plane as to the equatorial plane of Saturn. The x-axis of the Laplace place is taken as the node of 
the Laplace plane on the equatorial plane of Saturn. The obliquity amplitudes and resonance factors 
{f r j) °f tlie- solid Cassini state model presented in Section 3 are given in the last two columns. 
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1 


0.3197 


-0.00893124 


-703.51 


160.691 


0.1199 


1.38 
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0.0150 


-0.00192554 


-3263.07 


102.230 


0.0009 


1.06 
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0.0129 


0.42659824 


14.73 


292.867 


-0.0120 
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0.0022 


-0.21329912 


-29.46 


222.920 


-0.0026 


1.18 



as for the periodic variations about its equilibrium rotation rate, but with different assumptions 
concerning the timescales involved (long timescales for the obliquity and short timescales for the 
length-of day variations) . The influence of a global liquid layer on the LOD variations of Titan and 
on the librations of the Galilean satellites was investigated in [9] and pQ, respectively. Here, we 
extend their method, considering the appropriate timescales, to the Cassini state. By considering 
the gravitational and pressure torques between the different layers of the satellite, we show that the 
orientation of the rotation axis given in [6] can be partially reconciled with the moment of inertia 
given in [3]. 

2 The Cassini state for a solid Titan 

In the classical Cassini state, the rotation axis, the normal to the orbit, and the normal to the 
Laplace plane of a synchronous solid satellite remain in the same plane since the rotation axis has 
the same constant precession rate as the normal to the orbit about the normal to the Laplace plane, 
which is, by definition, the mean orbital plane of the satellite. The obliquity is then the constant 
angle between the rotation axis and the normal to the orbit. 

To be able to develop the Cassini state in the presence of a liquid subsurface ocean, we first 
present the generalization of the Cassini state of a solid Titan to a multi-frequency node precession, 
following [3J. Neglecting wobble, the rotation axis coincides with the principal axis of the polar 
moment of inertia C and the angular momentum equation is 

riTn{s A n), 

^MR 2 (-C 20 + 2C 22 )n T = ^(C-A)n T , (1) 

where s = (s x ,s y ,s z ) and h = (n x ,n y ,n z ) are the unit vectors along the rotation axis and the 
normal to the orbit, expressed in coordinates (x, y) of the Laplace plane and z along the normal. In 
addition, M and R are the mass and mean radius, C 2 q and C 22 are the second-degree gravity field 
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coefficients, A is the smallest moment of inertia, and nx is the mean motion, equal to the rotation 
rate, of Titan. The right-hand member is the gravitational torque (averaged over the orbital period) 
exerted by Saturn. The vectorial equation is then projected on the Laplace plane, where the motions 
of the projected rotation axis and orbit normal are easy to parametrize 

f (2) 

with S = Sx + Isy and N = n x + In y , I = \J — 1. Equation ^ is correct up to the first order in small 
obliquity, eccentricity and inclination (see Eq. (46) of [3] in which a sign has been corrected). We 
first assume that the orbital precession N, which causes the rotation axis precession, is zero in order 
to get the free spin precession mode whose frequency ujf depends only on the physical properties of 
the satellite 

= (3) 
To compute the forced solution, the orbital precession is written as a series expansion 

N = J^sini i e I to t4 ^-*/ 2 \ (4) 

3 

where the inclination amplitudes ij , frequencies Clj , and phases jj associated with the node preces- 
sion of the orbit with respect to the Laplace plane have been taken, up to j = 4, from [11] and are 
given in Table 1. The parameters Cli and i\ are the main precession rate (with a period of 703.51 
years) and the small mean inclination (0.3197°) of the orbital plane, considered in the classical 
Cassini state. The forced solution of Eq.Q is then 

S = ^sin(i i + £j V { ^ t+7 ^ /2) , (5) 

3 

where, correct up to the first order in Ej and ij, the obliquity amplitude £j associated with the 
frequency Clj is given by 

Sj = >j(lj . . (6) 

{u f + flj) 

The obliquity e at any time is the non-constant angle between the rotation axis and the normal to 
the orbit and oscillates between two extreme values £ m i n an d £ ma x such that 

e ^ sme=||5-JV|| (7) 
£ min = 2 max{|ej|} - V] \ej\ < e < V] \ej\ = e max . (8) 

3 3 

In this generalization, the normal to the Laplace plane, the normal to the orbit, and the rotation 
axis are not coplanar and the small angular deviation 5 of the spin axis with respect to the plane 
defined by the other axes is 

5 = sin<5 = (n x s y - n y s x )/\\N\\. (9) 
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Since the frequency ujf is positive, ujj + Clj can be close to zero if Clj < 0, and the coefficient Sj 
can be amplified by a resonance. For Clj < 0, we define "resonance factors", frj, which describe the 
amplification of Ej caused by the resonance and are close to one far from resonance, given by 

frj = max(|wj|, |f2j|)/|cjj + Clj\ (10) 

For C 20 = -31.808 x 10~ 6 , C 2 2 = 9.983 x 10~ 6 , and C = OMUMR 2 (less et al. 2010), the 
resonant factors frj are close to one because the free mode period is 191.91 years. Therefore, no 
significant resonant amplification (see Table 1) occurs and the obliquity variations are small. The 
obliquity e of about 0.12°, and the maximal deviation 5 of about 0.02° (Fig. 4) are inconsistent with 
the results of [6]. 



3 Titan with a liquid ocean 

3.1 A new Cassini state solution 

We now assume that Titan consists of four homogeneous layers : an ice shell (sh), a liquid ocean 
(o), an ice mantle (m), and an ice/rock core (c). The solid layer composed of the mantle and the 
core is also called interior (in) hereafter. The shell and the interior are considered to behave rigidly. 
For the solid layers, the angular momentum vector Hi can be written as the product of the polar 
moment of inertia C; and the rotation vector n^s;. The angular momentum equations take the form 



dsjh 
it 
dH n 



nrCsh—^r — f sh)ext + f p shext + f sh)int + T p shint , (11) 



dt 



dSi 



n TCi n — Ti n:ex t + ^ next + f i n ,int + ^ n ,int^ (13) 

where the external gravitational torque exerted by Saturn on layer (/), r^ ext , is equal to — f v (r A 
PjVW) dV, where W is the gravitational potential of Saturn averaged over the orbital period, r is 
the position vector of the points located inside the volume V\ of layer (/), and T^ int is the internal 
gravitational torque due to layers with a different orientation than layer (I). The liquid ocean induces 
pressure torques at the interfaces between the solid layers and the ocean. As we shall see later, the 
pressure torques can be interpreted as a modification of the external and internal gravitational 
torques and are therefore denoted by T p lext and rf int . 

The external gravitational torques on the shell and the interior are written as in Eq.Q for the 
solid case 

3 

F 8 h,ext = n T K sh {s sh A n), w\th.K sh = -(C sh - A sh )n T , (14) 

3 

T in ,ext = n T Ki n (s in A n), with K in = ^(Cin - A in )n T . (15) 

We divide the ocean into a top part and a bottom part, respectively, above and below an arbitrary 
chosen sphere inside the ocean. The pole axes of the top and bottom parts are those of the shell 
and the interior, respectively, and 

F ,ext = n T K 0jt (s s h A h) + n T Ko,b(sin A n), 

3 3 

K o,t = -{C 0>t - A 0:t )n T , K 0jb = -(C 0tb - A ofi )n T: (16) 
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3.1 A new Cassini state solution 



where (C ,t — Aj,t) and (C ,6 ~~ Aj,b) are the moment of inertia difference of the top and bottom part, 
respectively. For the long timescales considered here, it is a very good approximation to assume that 
the fluid is in hydrostatic equilibrium, therefore W induces a pressure P ex t in the ocean such that 
VP ex t = — p VW. The modification of the external torque on the interior because of the pressure 
is, applying Gauss' theorem, 

f?n,«* = I (rA Po VW)dV, (17) 
JV in 

where r is the position vector of the points located inside the volume Vi n of the interior. It then 
follows that 

?in,ext = n T K 0)b (s in /\n). (18) 

In the same way, the modification of the external torques on both the shell and the ocean because 
of the pressure effect are 

T P ,h,ext = n T K 0:t (s sh Ah), (19) 
f v oext = -n T K 0:b (s in An) - n T K 0)t (s sh An). (20) 



By using Eq. (16), we see that T 0fiX t + T p oext = 0, which is a consequence of hydrostatic equilibrium. 

Two layers with different obliquities and coincident axes of the moment of inertia B (on average, 
as a consequence of the Cassini state) exert a gravitational torque on each other tending to align 
their pole axis with each other. This internal gravitational torque on any layer (I) is given by 

Ti M t = ~ [ (?ApiV$)dV, (21) 
Jv t 

where is the internal gravitational potential, averaged over the orbital period, of the layers with 
a different obliquity than layer (/). From [8], we express the torque on the interior due to the shell 
and the top ocean and the torque on the shell caused by the misalignment with the bottom ocean 
and the interior as 

?in,int = -(8nG/5)(C in - A in ) [p s h(a sh - a + (3 sh /2 

-Po/2) + Po {a + /3 D /2)] (s sh A s in ), (22) 

f S)in t = (87rG/5) \{C in - A in ) + {C Qyb - A 0>b )] X 

[p sh (a sh -a + f3 sh /2 - f3 /2)](s sh A s in ), (23) 

where a and f3 are the polar and equatorial flattenings, defined as the relative differences ((a + 
b)/2 — c)/(a + b)/2) and (a — b)/a, respectively, with a > b > c the radii in the direction of the 
principal axes of the layers. For the ocean, we sum the torques on the top and bottom parts 

r 0) int = (8irG/5)(C in - A in )[p (a + f3 /2)](s sh /\s in ) 
-(8ttG/5)(C 0)6 - A, )6 ) x 

[psh{a sh -a Q + p sh /2 - /3 /2)](s sh A s in ). (24) 

The potential <3? induces a pressure Pi n t in the ocean (VPj n t = —p V(j)) that modifies the internal 
gravitational torque on the interior 

VLint = I (r/\PoV$)dV. (25) 
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Because of the similarity between Eq.([25j) and Eq.(21) and by using Eq.(23), we have 

rf„,int = -(8nG/5)(C 0>b -A 0ib )[ Psh (a sh -a + f3 sh /2 

-Po/2) + Po (a + Po/2))(§ sh A s in ). (26) 

In the same way as for the interior, we have 

? P sh ,int = (^G/5)[(C in -A in ) + (C , b -A , b )\ 

x [p {a + p /2)] (s sh A s in ) (27) 

for the shell. 

Therefore, T shjint + T p sh int + Ti n ,int + T p in ^ int = 0. Since the sum of all internal torques must be 

zero, we have that r 0j j nt + T p int = and that the total torque on the ocean is zero. Therefore, 
we only need to consider the angular momentum equations for the shell and the interior, given by 



Eqs.(ll) and (13). By introducing 

K 'sh = K sh + Ko,U (28) 
K in = K in + K o,fe) (29) 

K = -[(8vrG)/(5n T )][(Q n - A in ) + (C C)6 - A 0yb )\ x 

[ Psh (a sh -a Q + p sh /2 - p o /2) + p (a + p o /2)], (30) 

the system of angular momentum equation for the shell and the interior, projected onto the Laplace 
plane, is 

Csh^ = lK' sh (N-S sh )-IK(S in -S sh ), (31) 
CW^p = lK' in (N-S in ) + IK(S in -S sh ). (32) 

The two free modes of the system correspond to 'coupled' (with the frequency w+) and 'decou- 
pled' (w_) modes in which the shell and the interior oscillate in the same and opposite directions, 
respectively, such that 

w± = -(Z±VA)/(2C in C sh ), (33) 

Z = K{Ci n + C s h) — C s hK in — CinKgfr, 

A = -4C in C sh (-K(K in + K sh ) + K in K sh ) + Z 2 . 

For the orbital precession given by Eq.Q, we have, correct up to the first order in inclinations and 
obliquities, the forced solution for the spin positions 

Ssh = ^fe + ^.^e 7 ^'^-^, (34) 
j 

i 

£ _ ijQj(K(C s h + Cin) — C sn K in — CinC sn Qj) (36) 
C in C sh (uj + + (lj)(u- + (lj) 



ijQj(K(C s h + Ci n ) — Ci n K sh — Cin C s h £lj) 

C in C sh (uj + + + Clj) 



(37) 
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3.2 Numerical results 



Table 2 Size and density of the four internal layers of the Titan models consistent 
with the given constraints on mass and radius : The ocean thickness and the densities of the 
core are calculated for the given values of the other interior parameters. Within a given set denoted 
by curly brackets, the values are equally spaced. 



layer Thickness(km)/Radius(km) Density (kg m 3 ) 



ice shell 


h s h 


= {5}, {10, 20,. 


..,200} 


{800,900, ...,1200} 


ocean 


K = 


= 5-570 




{1000, 1100, ...,1400} 


ice mantle 


Rm 


= {2000,2025, 


..,2550} 


{1200, 1300, 1400} 


ice/rock core 


R c - 


= {1400,1410,. 


.,2200} 


2469 - 3176 



For comparisons with the observation, we define the shell obliquity (e s h), its minimum and maximum 
values (£max,sh,£min,sh), the deviation (5 s h), and two resonance factors (fff), which describes the 
amplification of £j :S h/in if &j < an d 0J± + is close to zero, in the same way as for the solid case 
(Eqs. Q, (§, © and ph). 



3.2 Numerical results 

The IAU orbit orientation used in [G] is less precise than their determination of the rotation axis 
orientation and is the main source of error in the obliquity and deviation calculation. By comparing 
the IAU orbit orientation with those derived from the ephemerides of the JPL HORIZONS system 
and of the TASS1.6 ([H]), we estimated the obliquity to be e = 0.32 ± 0.02° and the deviation to 
be <5 = 0.12 ±0.02°. 

We evaluate the solution for a representative range of hydrostatic interior structure models with 
homogeneous layers of constant density that are constrained by the observed mass and radius (see 
Table 2). The moments of inertia of the layers depend on the polar and equatorial flattenings due 
to rotation and static tides that are computed from the Clairaut equation, assuming hydrostatic 
equilibrium. Among these models, we retain those that have a moment of inertia consistent at the 
3cj level with the estimated value of 0.3141 ± 0.0005 of @]. 

The periods of the free modes (T± = 2ir/u±) and the shell obliquity amplitudes £j, s h as a function 
of the moment of inertia for the interior models are shown in Figs. 1 and 2, respectively. Far from 
the resonances (/r - i i 24 — 1), the shell obliquity is below 0.15°. Therefore, a significant resonant 
amplification of at least one of the Ej S h is needed to obtain a time- variable obliquity that can be as 
large as the observed one ([£ m ax,sh'i £min,sh] H [0.32° — 3<r; 0.32° + 3a] == 0). Some interior models have 
T + close to the opposite of the period of Cl\ or T_ very close to the opposite of the period of and 
have a resonant amplification of £\ iS h or £4^, respectively, that enables a sufficiently large obliquity 
value. These interior models are indicated by red and green dots in Figs. 1 and 2, which show that 
a specific interior model cannot be close to both resonances. Given our range of interior models (see 
Table 2), some models have the same moment of inertia, but not necessarily the same free mode 
period, and form vertical lines in Fig. 1. Because of the large number of retained interior models, 
these vertical lines cannot be easily distinguished, except for those appearing as peaks on the graph 
of T + and corresponding to models with the same h s , R c , p m , p , and p s , such as p Q = p m , and with 
different R m . One can show that T + is approximately proportional to for these models and is 
larger for larger R m . The interior models at the top of the peaks (red dots), hereafter called class 1 
models, are of interest because their shell obliquity can reach the observed value of the obliquity, as 
a result of a resonant amplification of £\. s h- These models are characterized by a very large interior 
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Figure 1 - Periods of the 'coupled' (T + ) and 'decoupled' (T_) free modes as a function of the 
moment of inertia, for our range of interior models. The blue horizontal lines represent the opposite 
of the first and the fourth orbital periods. Some interior models have T + close to 703.51 years and 
other interior models have a period T_ very close to 29.46 years, leading to a resonance. The red 
and green dots are the interior models that, thanks to a resonant amplification of E\ B or e^ s , can 
reach a shell obliquity as large as the observed one (see also Fig. 2). 



radius of 2525-2550 km, a thin ice shell (5-30 km thickness), and a thin ocean (5-45 km thickness). 
Interior models that can reach the observed obliquity through a resonant amplification of e^sh (class 
2 models) are more diverse than the models of the first class. They are characterized by a thicker 
ice shell (from 150 to 200 km), an ocean of from 5 to 425 km in thickness, an ice mantle of from 10 
to 590 km in thickness and a core radius of between 1800 and 2070 km. Table 3 presents two such 
resonant models. Model 1 is in resonance with the main period of precession (frf = 5.31 > 1) and 
model 2 is in a close resonance with the fourth period of precession (/rj = 171.11 S> 1). 

In addition to the obliquity, the deviation of the shell rotation axis from the plane defined by 
the orbit and Laplace plane normals (see Eq.Q) has to be consistent with the observations. With 
a small deviation of 0.03° at best, the class 1 models cannot explain the observed deviation of 
about 0.12°. The deviation reported by [6] can be seen as an averaged position of the rotation axis 
during the period of coverage of the analysed radar observations (from Oct 26 2004 to Feb 22 2007), 
during which model 2 has a deviation smaller than 0.12°. This is true for all the models of class 2 
with a positive £4^, because the deviation depends mainly on the ephemerides, particularly on the 
phases of the orbital precession. We would need to add about 15° and 45° to 74 to obtain the correct 
deviation at time Oct 26 2004 and Feb 22 2007, respectively. However, a comparison of TASS1.6 with 
a series expansion fitted on the ephemerides of the JPL HORIZONS system convinced us that the 
phase 74 in TASS1.6 is correct up to a few degrees. The series expansion given in [11] is incomplete 
and additional terms with smaller amplitudes exist (up to j = 21, [10J). With those additional 
terms, some internal models that are consistent with the moment of inertia can account for both 
the estimated obliquity and deviation, thanks to a combination of a very close resonance between 
ui + and £7i4 = — 27r/578.73y, besides the further resonance with (1%, discussed before. However, this 
very close resonance seems implausible because of the small fraction of suitable internal structure 
models. 
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3.2 Numerical results 



K 1,sh ^2,811 "S.sh l £ 4,5hl 

io ] 

10" 2 

io' 4 
io" G 

0.34 0.342 0.34 0.342 0.34 0.342 0.34 0.342 
C/MR 2 C/MR 2 C/MR 2 C/MR 2 

Figure 2 - Shell obliquity amplitudes £j jS h as a function of the moment of inertia, for our range 
of interior models. £\ iS h and E2, s h ar e positive whereas is negative because of the positive sign 
of Q3. e^^h can be either negative or positive since (w_ + O4) can be positive or negative. The 
red/green dots have the same meaning as in Fig. 1. 




Table 3 - Internal structure parameters, obliquities, free modes periods, and resonant 
factors for two resonant models : Models 1 and 2 are in resonance with the first and the fourth 
period of precession, respectively. 







Model 1 


Model 2 


hshi hoi 




10, 15, 2550, 1680 


170, 105, 2300, 1890 


Psh, Po, 


Pmi Pc 


1200, 1400, 1400, 3142 


1100, 1400, 1400, 2758 


£l,sh, 


£l,in 


0.3198, 1.3816 


0.0999, 0.2528 


£2,sh, 


£2. in 


0.0007, 0.0032 


0.0007, 0.0016 


£3,sh, 


£3, in 


-0.0043, -0.0126 


-0.0102, -0.0125 




£&,in 


-0.0004, -0.0023 


0.2300, -0.0270 


T+, 


r_ 


572.019, 2.342 


300.391, 29.285 



fr[ w , fr™ 9 , fr" g 5.31, 1.21, 1.05 1.75, 1.10, 1.11 

frf ec , M ec , frf° 1.00, 1.00, 1.09 1.04, 1.01, 171.11 
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v 




x 



Figure 3 - Projection of the normal to the orbit (blue) and of the rotation axis on the Laplace 
plane (black for the solid case, red for model 1 and green for model 2) over the period of the main 
precession, beginning at J2000 (thin curves) and over the observation period of [6] (thick curves). 
The "+" marker is the projection of the rotation axis measured by [6]. (Unit of the graph is radian.) 




S 10 IS 20 25 30 ' 5 10 J5 20 25 30 

Time (years) Time- (years) 



Figure 4 - Shell obliquity e s h (left) and deviation 6 s h (right) over 30 years beginning on J2000 
(black for the solid case, red for model 1 and green for model 2). The vertical blue lines show the 
observation period of [6] and the blue boxes represent our estimated obliquity (0.32° ± 3<r) and 
deviation (0.12° ±3<r). 
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4 Discussion and perspectives 

Since we have found a better consistency between the measured and computed positions of the 
rotation axis with our model than in the solid case, assuming that Titan is locked in the Cassini 
state, our study is an indication of a possible water ocean layer beneath the surface of Titan. 

For two classes of interior models with a liquid water ocean beneath an ice shell, the obliquity 
computed with our new Cassini state model can be as large as the observed one. However, we 
rejected the first class of models with a very thin ice shell and a shallow liquid ocean since a thin 
shell implies high temperatures and therefore, according to the water phase diagram, a thick ocean. 
On the other hand, models of class 2 are possible but can be quite different so that the estimated 
obliquity and moment of inertia do not provide accurate constraints on the interior structure. 

The theoretically predicted deviation of the rotation axis with respect to the plane defined by 
the orbit and Laplace plane normals is inconsistent with the observations. We have shown that the 
use of ephemerides extended to additional frequencies could solve this problem, but only for a small 
fraction of the rejected class 1 models. Another explanation is that Titan may be slightly offset from 
the Cassini state because of a recent excitation. However, other explanations extending our Cassini 
state model might be tested. We have neglected the viscosity at the ocean boundaries because we 
have assumed that the timescale on which the viscosity is effective is greater than the timescale of 
the node precession of Titan. We have also neglected the polar motion induced by the atmosphere 
and the resulting variations in the rotation axis orientation in space since they are expected to have 
a very small amplitude compared to Titan's obliquity. Although the deformation of the ice shell has 
not been included, our new Cassini state model is a first step towards a more realistic model of the 
obliquity of an icy satellite with an ocean. 
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